#include <stdio.h>
#include <math.h>

#include "timer.h"

#include <omp.h>

#define CHUNK 512

void ReadInputData(char* inputDataFileName, int* nx, int* nt)
{
	FILE* file;
	file = fopen(inputDataFileName, "r");

	if (file != NULL)
	{
		fscanf(file, "%d %d", nx, nt);
		fclose(file);

		printf("%d %d\n", *nx, *nt);
		return;
	}

	printf("Error openimg file\n");
}

void WriteOutputData(char* outputDataFileName, double data)
{
	printf("%f\n", data);

	FILE* file;
	file = fopen(outputDataFileName, "w");

	if (file != NULL)
	{
		fprintf(file, "%f", data);
		fclose(file);

		return;
	}

	printf("Error openimg file\n");
}

double GetMaxDeriation(double* expected, double* counted, int size)
{
	double max = 0;
	for (int i = 0; i < size; i++)
	{
		if (max < fabs(expected[i] - counted[i]))
		{
			max = fabs(expected[i] - counted[i]);
		}
	}
	return max;
}

double AccurateSolve(double x, double t)
{
	return x * x * x + t * t;
}

double f(double x, double t)
{
	return 2 * t - 6 * x;
}

int GetGreaterPow(double x)
{
	return pow(2, floor(log(x) / log(2.0)) + 1);
}

void Initialize(double* &a, double* &b, double* &c, double* &rightPart, double* &x, int nx, int nt, double* &aReduce, double* &bReduce, double* &cReduce, double* &rightPartReduce)
{
	int n = GetGreaterPow(nx);
	int reduceStepCount = log((double)n) / log(2.0);

	a = new double[n - 1];
	b = new double[n - 1];
	c = new double[n - 1];

	rightPart = new double[n - 1];
	x = new double[n + 1];

	int sum = 0;
	for (int k = 1; n > 1; n /= 2, k++)
	{
		sum += (n - 1);
	}

	aReduce = new double[sum];
	cReduce = new double[sum];
	bReduce = new double[sum];
	rightPartReduce = new double[sum];
}

void GenerateLeftPart(double* a, double* b, double* c, int nx, int nt)
{
	register double acCoef = -nx * nx;
	register double bCoef = 2 * nx * nx + nt;
	int i;

#pragma omp parallel shared(a, b, c) private(i)
	{
#pragma omp for schedule(static, CHUNK) nowait
		for (i = 0; i < nx - 1; i++)
		{
			a[i] = acCoef;
			b[i] = bCoef;
			c[i] = acCoef;
		}
	}

	int greaterPow = GetGreaterPow(nx);

	if (greaterPow > nx)
	{
		for (int i = nx - 1; i < greaterPow - 1; i++)
		{
			a[i] = 0.0;
			b[i] = 1.0;
			c[i] = 0.0;
		}
	}
}

void GenerateRightPart(double* rightPart, double* lastGrid, int nx, int nt, int t)
{
	int i;

	double time = (double)t / (double)nt;

#pragma omp parallel shared(rightPart, lastGrid) private(i)
	{
#pragma omp for schedule(static) nowait
		for (i = 0; i < nx - 1; i++)
		{
			rightPart[i] = lastGrid[i + 1] * nt + f((double)(i + 1) / (double)nx, time);
		}
	}

	lastGrid[0] = AccurateSolve(0.0, time);
	lastGrid[nx] = AccurateSolve(1.0, time);

	int greaterPow = GetGreaterPow(nx);

	if (greaterPow > nx)
	{
		lastGrid[greaterPow] = 0; 
		rightPart[nx - 1] = AccurateSolve(1., time);

		for (i = nx; i < greaterPow - 1; i++)
		{
			rightPart[i] = 0.0;
		}
	}
}

void Terminate(double* a, double* b, double* c, double* rightPart, double* x, double nx, double* aReduce, double* bReduce, double* cReduce, double* rightPartReduce)
{
	delete[] a;
	delete[] b;
	delete[] c;

	delete[] rightPart;
	delete[] x;

	delete[] aReduce;
	delete[] bReduce;
	delete[] cReduce;
	delete[] rightPartReduce;

}

void CircularReduce(double* a, double* b, double* c, double* rightPart, double* x, int size, double* aReduce, double* bReduce, double* cReduce, double* rightPartReduce)
{
	register int i, i2, id;
	int threadNum = 16;

	register double alpha, betta;

	int n = size;
	int reduceStepCount = log((double)n) / log(2.0);

#pragma omp parallel shared(aReduce, bReduce, cReduce, rightPartReduce, a, b, c) private(i)
	{
#pragma omp for schedule(static) nowait
		for (i = 0; i < n - 1; i++)
		{
			aReduce[i] = a[i];
			cReduce[i] = c[i];
			bReduce[i] = b[i];
			rightPartReduce[i] = rightPart[i];
		}
	}

	register int stringNumberK_1 = 0;
	n /= 2;

	for (int k = 1; n > 1; n /= 2, k++)
	{
#pragma omp parallel shared(aReduce, bReduce, cReduce, rightPartReduce, k, stringNumberK_1, n) private(i, i2)
		{
#pragma omp for schedule(static)
			for (i = 0; i < n - 1; i++)
			{
				i2 = 2 * i;

				alpha = aReduce[stringNumberK_1 + i2 + 1] / bReduce[stringNumberK_1 + i2];
				betta = cReduce[stringNumberK_1 + i2 + 1] / bReduce[stringNumberK_1 + i2 + 2];

				aReduce[stringNumberK_1 + (2 * n - 1) + i] = -aReduce[stringNumberK_1 + i2] * alpha;
				bReduce[stringNumberK_1 + (2 * n - 1) + i] = bReduce[stringNumberK_1 + i2 + 1] - cReduce[stringNumberK_1 + i2] * alpha - aReduce[stringNumberK_1 + i2 + 2] * betta;
				cReduce[stringNumberK_1 + (2 * n - 1) + i] = - cReduce[stringNumberK_1 + i2 + 2] * betta;
				rightPartReduce[stringNumberK_1 + (2 * n - 1) + i] = rightPartReduce[stringNumberK_1 + i2 + 1] - alpha * rightPartReduce[stringNumberK_1 + i2] - betta * rightPartReduce[stringNumberK_1 + i2 + 2];
			}
		}
		stringNumberK_1 += (2 * n - 1);
	}

	int start = size / 2;
	for (int k = reduceStepCount; k > 0; k--, n *= 2, start /= 2)
	{
		for (int i = 0; i < n; i++)
		{
			i2 = i * 2;
			id = start + i2 * start;

			x[id] = (rightPartReduce[stringNumberK_1 + i2] - aReduce[stringNumberK_1 + i2] * x[id - start] - cReduce[stringNumberK_1 + i2] * x[id + start]) / bReduce[stringNumberK_1 + i2];
		}

		stringNumberK_1 -= (n * 4 - 1);
	}
}

void solve(double* a, double* b, double* c, double* rightPart, double* x, int nx, int nt, double* aReduce, double* bReduce, double* cReduce, double* rightPartReduce)
{
	register int i;
	register int greaterPower = GetGreaterPow(nx);

	//if(greaterPower < 17000) omp_set_num_threads(2);
#pragma omp parallel shared(x) private(i)
	{
#pragma omp for schedule(static) nowait
		for (i = 0; i < nx + 1; i++)
		{
			x[i] = AccurateSolve((double)i / (double)nx, 0);
		}
	}

	GenerateLeftPart(a, b, c, nx, nt);

	for (int j = 1; j < nt + 1; j++)
	{
		GenerateRightPart(rightPart, x, nx, nt, j);

		CircularReduce(a, b, c, rightPart, x, greaterPower, aReduce, bReduce, cReduce, rightPartReduce);
	}
}

int main(int argc, char** argv)
{
	int nx, nt;
	double* a = NULL, *b = NULL, *c = NULL;
	double* rightPart = NULL;
	double* expected; 
	double* calculatedSolve = NULL;
	double maxDeviation;

	double* aReduce = NULL; 
	double* bReduce = NULL;
	double* cReduce = NULL;
	double* rightPartReduce = NULL;

	char* inputFileName;
	char* outputFileName;
	char* spentTimeFileName;

	if (argc != 4)
	{
		printf("Wrong parameters count\n");
		return 0;
	}

	inputFileName = argv[1];
	outputFileName = argv[2];
	spentTimeFileName = argv[3];

	ReadInputData(inputFileName, &nx, &nt);

	Initialize(a, b, c, rightPart, calculatedSolve, nx, nt, aReduce, bReduce, cReduce, rightPartReduce);

	Timer timer;
	timer.start();
	solve(a, b, c, rightPart, calculatedSolve, nx, nt, aReduce, bReduce, cReduce, rightPartReduce);
	timer.stop();

	expected = new double[nx + 1];
	for (int i = 0; i < nx+1; i++)
	{
		expected[i] = AccurateSolve((i) * 1.0 / (double)nx, 1.0);
	}

	maxDeviation = GetMaxDeriation(expected, calculatedSolve, nx+1);

	delete[] expected;
	Terminate(a, b, c, rightPart, calculatedSolve, nx, aReduce, bReduce, cReduce, rightPartReduce);

	WriteOutputData(outputFileName, maxDeviation);
	WriteOutputData(spentTimeFileName, timer.getElapsed());

	return 0;
}